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ABSTRACT 


In this thesis work some applications of Wavelet transform in the area of Image Recon- 
struction from Projections have been studied. The process of data collection in computer 
aided tomography, introduces noise in the projections. One of the major sources of noise is 
the statistical nature of photon emission and detection. The wavelet transform is used to 
reduce the effect of noise in the reconstructed image by preprocessing the noisy projections. 
This method is compared w r ith the existing technique and its advantages and disadvantages 
are discussed. It is shown that reconstruction of images from the ID wavelet transform 
of the projections gives the 2D wavelet transform of the original image. The 2D wavelet 
corresponding to this wavelet transform is shown to be the inverse Radon transform of 
the ID wavelet. These two results are combined to produce the noise suppressed analysis 
images from the noisy projections. The above results are verified through simulations using 
computer generated images. 
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Chapter 1 
Introduction 


1.1 Review 

The problem of determining some property of the internal structure of some object without 
macroscopically damaging the object arises in diverse area such as medicine, electron mi- 
croscopy, stress analysis etc. Radon transform is the unifying mathematical framework for 
all such problems. In medicine, X-ray transmission tomography, also known as computer 
tomography (CT), is one of the important techniques used to obtain the x-ray attenuation 
coefficient distribution over a cross section of some portion of the human body from the 
projection profile that are measured at various angles. The process of' collecting projection 
data has many practical limitations, thereby, introducing error in the projection measure- 
ments. One such basic limitation to the accuracy of measurements in CT is the statistical 
nature of the photon production, its interaction with the matter and photon detection [3]. 
It is shown that this error can be modeled as Additive Gaussian noise with zero mean [13]. 
At high frequency, this noise dominates the signal. The reconstruction process, which in- 
volves filtering the projection with ramp filter, amplifies the high frequency content of the 
input signal. This effectively increases the noise in the reconstructed image. In practice 
this problem is tackled by using a lowpass filter in conjunction with the ramp filter. This 
low pass filter will have gentle roll off characteristic to avoid ringing. The problem with 
this spatially invariant filtering is, it suppresses the high frequency content of the signal 
along with the noise. Thus the improvement in noise performance comes at the expense of 
resolution. 

Recent lyfcfew methods using spatially variant filtering have been proposed in the lit- 
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erature [9, 10]. In [9], Sahiner et.al., propose a spatially varying filtering using time fre- 
quency distribution. Window Fourier transform is used to represent the signal in the 
time-frequency plane. In [10], the same authors propose a method using wavelets, where 
prior knowledge of the image is used to set some of the coefficients of the 2D wavelet 
transform of the image to zero and reconstruct the image imposing these constraints. 

One major problem in window Fourier transform is in choosing an appropriate window 
(the analysis filter). Too short a window may result in poor frequency resolution and 
too long may result in poor spatial resolution. This problem does not arise if we choose 
wavelet transform instead of window Fourier transform, because in wavelet transform, the 
window size varies as a function of scale. Also fast algorithms exist in computing the 
wavelet transform. With this as motivation we propose a method using spatially varying 
filtering in time-frequency plane using wavelet transform. Unlike in [10], we assume no 
prior knowledge of the original image. 

Wavelet transform of a given signal can be expressed as convolution of the given signal 
with rp a (—s), the reversed wavelet at a given scale ‘a’ , where ip(s) is an admissible one 
dimensional wavelet[5]. From the convolution property of the Radon transform, we know 
that the Radon transform of the two dimensional convolution of two signals is the same as 
the one dimensional convolution of the Radon transform of the two signals on the variable 
‘s’ [1]. Using this property it can be shown that the reconstructed images from the wavelet 
transform of its projections can be said to be two dimensional convolution of the original 
2D image with a 2 dimensional function 'F a (— x, — y) which is the inverse Radon Transform 
of the one dimensional analysing wavelet (rather its reflected and scaled version, to be 
precise). We show in this work that the two dimensional function $„(— x, — y) is a reversed 
and scaled version of a two dimensional function i.e., \k a (— x, — y) = 'F(— x/a, — y/a). We 
also show that ^(ar,y) is an admissible two dimensional wavelet [12]. 

1.2 Organization of Thesis 

• In chapter 2, we review the relevant results in Radon transform, Wavelet transform 
and Detection theory. 



Chp.l Introduction 


3 


• In chapter 3, we propose a method to suppress noise in the reconstructed image from 
noisy projections using wavelet transform. 

• In chapter 4, we prove that ^(x, y) is an admissible wavelet and hence show that the 
2 dimensional convolution is actually the 2 dimensional wavelet transform equation. 

• We test the above propositions on computer simulated projections of the synthetic 
images and present the results in chapter 5. 

• Suggestion for future work is presented in chapter 6. 

• We summarise our results and conclude in chapter 7. 

In this thesis work we use 71 to represent the set of real numbers and 2 to represent 
the set of integers. 



Chapter 2 

Mathematical Preliminaries 


In this chapter a very brief review of the results that will be used in this work are given. 
We direct the interested reader to texts and journal articles, for a detailed treatment on 
these topics. These references are listed in Bibliography. In section 1 we review Radon 
transform and some of its properties. In section 2, we present the relevant results in wavelet 
transform and multiresolution analysis. In section 3, we give two results from detection 
theory which will be used in this work. 

2.1 Radon Transform 

Image reconstructions from projections is the process of producing an image of a two 
dimensional distribution (of some physical property) from the estimates of its line integral 
along a finite number of lines of known locations. The basic mathematical framework 
common to a large class of such reconstruction problem was first developed by Johann 
Radon as early as in 1917. His work which has come to be known as Radon transform on 
Euclidean space (7£ n ), relates a function defined on n-dimensional Euclidean space 7Z n and 
its integrals over all the hyperplanes of dimension n-1. Thus if f is defined over a plane 
H 2 , like in 2 dimensional medical imaging, its Radon transform is determined by its line 
integrals over all lines passing through the function support. Knowledge of all possible 
line integrals constitute full knowledge of the Radon transform. In Practice we have only 
a finite number of samples of the Radon transform. The image reconstruction techniques 
deal with the problem of reconstructing the image from this partial information. 
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Figure 2.1: The coordinates in space domain and the line L representing the line of inte- 
gration 

Definition : Let (x, y ) designate coordinates of a point in a plane. Consider an arbitrary 
function /(x, y ) defined over some domain D in R 2 . Then its Radon transform (also known 
as the projections), g(s,0 ) is given by 

/ oo roo 

/ f(x,y)6(xcos0 -f ysinO — s) dx dy (2-1) 

-OO J — OO 

If f is continuous and has compact support, then g(s,0) uniquely determines /(x,y) [1]. 

In equation 2.1, d(.) is a Dirac Delta function. Thus for a given 6 , the equation 2.1 
is a line integral along the line xcos# + ysin# = s (figure 2.1). Associated with Radon 
transform is the Back projection operator which is given by 

b(x,y) = f g(x cos 6 -f y sin 0, 6) d0 (2.2) 

Jo 

Back projection represents the accumulation of the ray sums of all the rays that pass 
through the point (x,y). Backprojection operator is not the inverse of the Radon transform 
but the inverse Radon transform operation involves this operator. The inverse Radon 
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transform is, 

/(*,?) = 4^ [ JZ U I g{U)e 3i(xco,e ^ axne) did9 (2.3) 

where g{£,0) is the one dimensional Fourier transform of g(s, 9) with respect to ‘s’. We 
note that in equation 2.3, the inner integral represents the inverse Fourier transform of 
|£li(£>0)) i- e -> if represents the filtering of the projections g(s,9 ) with the ramp filter p(s) 
whose frequency response is |£|. The outer integral represents the backprojection operation 
on these filtered projection. 

As mentioned in chapter 1, the ramp filter amplifies the high frequency components of 
the projections. Since the noise strength in the high frequency region is generally more 
than the signal strength, it is common practice to use a low pass filter w(£) in conjunction 
with the ramp filter, to improve signal to noise ratio. Hence the filter |£| in equation 2.3 is 
replaced by h(£) = |£|u>(£). The low pass filter usually has a gentle roll off characteristics 
to avoid ringing at the edges. One such filter A(£), which exhibits good noise suppressing 
property is [2, 13], 

MO = l£l (« + (1 - <*) cos 2x£<£] for |£| < 1 /24 (2-4) 

where d is the sampling interval along the variable ‘s’. This filter is called generalized 
Hamming filter. But this space invariant lowpass filtering results in loss of resolution, 
since the filter h{ 0 smoothes out the edges in the image. We will come back to this 
problem in chapter 3. 

Some properties of the Radon transform which we will refer to in this work are given be- 
low. In the following properties g(s,9), gi(s,0) and g 2 (s,0) represent the Radon transform 
of the functions f{x,y), fi(x,y) and f 2 (x,y) respectively. g(£,0) is the Fourier transform 
of g(s,9 ) and is the 2 dimensional Fourier transform of f(x,y). 

1. Projection Slice theorem : The one dimensional Fourier transform on s of g(s,9) is 
equal to the central slice, at an angle 9 of the 2 dimensional Fourier transform of the 
object f(x,y). 

g((,9) = F(£ cos0,£sin$) 

2. Radon transform is linear, i.e., given the scalars, ai,a 2 € ft, we have, 

R T 

aifi(x,y) F a 2 f 2 (x,y) *—> aigi(s,6) + a 2 g 2 (s,0) 


(2.5) 
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3. Convolution property : The Radon transform of the 2 dimensional convolution of 
two 2-D signals f\{x,y) and f 2 (x,y) is the 1 dimensional convolution of the Radon 
transform of the two signals. 

/i( z >y) ®2 Mx,y) 9i(s,0) ®1 g 2 (s,e) (2.6) 

where <S>i represents one dimensional convolution on ‘s’, ®2 represents two dimen- 
sional convolution and R.T represents the Radon transform operation. 

In Computer Tomography, the linear attenuation coefficient, fi(x ,y), of the material to 
x-ray is the quantity we are trying to estimate. We use the data collected through x-rays 
which gives an estimate of the line integrals. 



where I 0 is the input intensity (number of photons per sec per unit area), I is the observed 
intensity after the beam passes through the material and du is the increment length along 
the line L (figure 2.1). In actual measurement due to practical limitations such as the 
statistical nature of the photon emission and detection, beam hardening etc.[3] noise is 
added to the projections. In [13] Shepp et al., show that noise introduced due to the 
statistical nature of photon can be modelled as an additive Gaussian noise with zero mean. 

Proofs of the above properties, are given in [1, 2]. The computer implementation 
of discretized reconstruction algorithm is discussed in detail by Rowland in [4]. In [3], 
Herman presents a detailed discussion of various reconstruction techniques as related to 
computer tomography and its practical limitations. For an overview of Radon transform, 
its properties and an English translation of Radon’s original paper, refer [1]. 

2.2 Wavelet Transform 

Wavelets are relatively recent development in applied mathematics. The concept of wavelets 
can be viewed as a synthesis of ideas which originated in diverse areas such as, physics, 
engineering, and pure mathematics. As a result, wavelets appeal to engineers as well as to 
scientists. Image coding, speech processing, coherent state study in physics, are but few of 
the applications of wavelets. 
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The Fourier transform of a given signal gives a frequency content of the signal. But it 
fails to provide information about the time evolution of the frequency content, i.e., how the 
frequency contents of the given signal change with time. This time-frequency localization 
is achieved by wavelet transform. The wavelet transform of signal evolving in time depends 
on two variables, scale (a) and time (b). 

= (2 ' 7) 
where ip(t) is the mother wavelet which satisfies the admissibility condition given by 

C V = 2tt / + °° di < oo (2.8) 

y-oo |f| 

where ip{£) is the Fourier transform of This requires %l){t)dt = 0. When the 

wavelet satisfies the above condition, the inverse wavelet transform exists and is given by 

m = k £ £ ("vx* *> * (¥)& db (2 - 9) 

This scaling in the wavelet transform helps to zoom into the short lived high frequency 
signal or zoom out to the lengthy low frequency signal, by varying the width of the support 
of the function t/>(£) [5, 6]. This property of the wavelet transform distinguishes it from 
other time frequency analysis techniques such as window Fourier transform. The discrete 
wavelet transform is obtained by sampling the parameters a and b in the following manner, 
a = a™ and b = nb 0 a™ with a 0 > 1 and 6o > 0 

In general the continuous and discrete wavelet transform are redundant representation 
of a given signal. For some very special choice of tl>,a 0 and 6 0 , the redundancy in the 
wavelet representation is removed. It has been proved that if we choose % = 2 and 6q = 1 
then {tfimn, : m,n € Z} constitute an orthonormal basis of L 7 (1Z) ([ 5 ]). This removes the 
redundancy in the representation. In this case we define V’mn = 2' -m/20(2 _ n ) The 
wavelet transform pair is now given by, 

P-10) 

and 

fit) = Yj /m»Vw(0 

m,n 


( 2 . 11 ) 
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It has been shown by Mallat[l4] that orthonormal wavelets are closely related to Mul- 
tiresolution Analysis (MRA). A sequence of vector spaces {V : : j € Z} is said to be a 
Multiresolution Analysis of L 2 (7Z) if they satisfy the following properties and the theorems 
stated below. 

i. • . • C V 2 C Vi C Vo C K_i C VL 2 C . . . 

ii. U je2 Vj = l\ n) 
i”- fW V, = (0} 

iv. if /( t) € Vj then f(2H) E V 0 

v. if f(t) e V 0 then f(t - k) e V 0 for k G Z 

Theorem : 1 For a given multiresolution analysis of L 2 {7i), there exists a function 
such that { <f>mn(t ) : n € Z) forms a basis of V m , where <f> mn {t) = — n). The 

function 4>{t) is unique and is the called scaling function. 

Theorem : 2 Associated with each MRA, is a unique orthonormal wavelet ip(t). 

MRA, besides providing a framework for understanding of wavelet bases, has led to 
the development of fast algorithms to compute Discrete wavelet transforms and Inverse 
discrete wavelet transforms by indicating the relation between the wavelet transform and 
subband techniques[5]. This relation between discrete wavelet transform(DWT) and the 
subband technique also establishes another important property of the DWT, namely, the 
DWT of a signal divides the frequency spectrum of the given signal into frequency bands 
with constant relative bandwidth. More specifically, the DWT at a scale 2 m , gives the 
frequency content of the signal in the range of < f < 2 ^r> T being the sampling 

interval. Hence the energy of the DWT coefficients at a scale 2 m gives an estimate of the 
signal energy in that frequency band. 

Figures 2.2 and 2.3 show the subband implementation scheme to compute the DWT 
and the IDWT . In the figure 2.2, h(n) = h(-n) and <?(n) = g(-n), where h(n) and g(n) 
are respectively the discrete low pass and the high pass filters associated with i/> ( < ) and 
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Figure 2.2: Computing wavelet coefficients through subband technique 
4>{t) given by equations 2.12 and 2.13. 


H n ) = 'fwmt-n)* 

(2.12) 

1 

■to 

■to 

II 

^5) 

(2.13) 


for n e 2 [5, 14, 15]. 

The sequence co n , in the figure 2.2 and 2.3 is the inner product of /(f) with 4>On{t) at 
the fine resolution level 2°. {d mn : m = 1,2..., M] represent the wavelet coefficients 
which are the detailed signals at the given resolution level 2~ m (scale 2 m ). They represent 
the difference in information between the blurred signals C( m _ X ) n at resolution 2~h 7l_1 ) and 
c mn at resolution 2~ m . Hence they are also called as difference signals. Similarly c mn 
represent the blurred signal at the given resolution level 2 -m , which represent the blurred 
version of the original discrete signal con- 


Mu, 



Figure 2.3: Inverse wavelet transform through subband techniques 











Chp.2 Mathematical Preliminaries 


11 


2.3 Detection Theory 

Consider the following detection problem, 

y — x + w : Hi 
y = w : HO 

where w is the standard normal random variable N( 0 ,<7 2 ) and £ is a deterministic signal 

with values in the range [R1,R2], Rl < R2. The problem of determining whether y is due 

to signal (hypothesis HI) or due only to noise (hypothesis HO) is optimally solved through 

the likelihood ratio test when Rl,R2 are both positive[7j. In this case this is an Uniformly 

Most Powerful test (UMP). The likelihood ratio test is given by, 

hi 

> 

y < 7 

HO 

where 7 is the threshold. The threshold can be fixed, for e.g., by specifying the desired 
probability of false alarm. 

If 721 < 0 and R2 > 0, it can be shown that the UMP doesnot exist ([7]). Then we 
resort to some suboptimal tests such as the generalized likelihood ratio test (GLRT), given 

by 

Hi 

Ivl < t (2-14) 

HO 

The equation 2.14 will be used in our work. 



Chapter 3 

Wavelet Preprocessing of Noisy 
Projections 


In this chapter we discuss the problem of noise suppression and propose a possible solution 
through time variant filtering, within the following framework. 

• The noise is Additive White Gaussian. 

• Signals considered are spacelimited and resolution limited . Essentially bandlimited 
signals which are represented by finite number of samples belong to this class of 
signals [10, 11]. 

• No a priori knowledge about the nature of input image and its projections is assumed. 

3.1 Time variant filtering 

The fine scale wavelet transform coefficients (those corresponding to smaller m in 2 m ) rep- 
resent the localized high frequency contents of the signal. Hence, setting all the fine scale 
coefficients to zero will amount to low pass filtering the original signal. This process, which 
is equivalent to time invariant filtering, reduces noise, at the expense of edge information 
present in those high frequency components. If we choose to set only those wavelet co- 
efficients whose absolute values are less than a particular threshold to zero, then we can 
still preserve the edges while reducing the effect of noise [9] . Such selective zeroing of the 
wavelet coefficients, which is equivalent to time variant filtering, is justified by the GLRT, 
discussed in section 2.3. 
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Let y{s x ,0j),g(s„0 } ) and w(s,,9j) represent the discrete forms of the noisy projections 
noiseless projections [#(.s, 0)] and the white noise [ie(s, 0)). Here, s, € N = 
{0, D , . . . , (N — !)£>} and 9 : € 0 = {0, n/K, 2 -kJK, . . . , (K — 1 )7r//T}, where N and K 
are the number of samples in ‘s’ and 6 respectively and D is the sampling interval along 
‘s’. Hence, the discrete noisy projections can be written as, y(s,,9j) = g(s t ,9j) -J- w(s t ,0j) 
The problem is to estimate g(.) from y(.). We can represent the signal y(.) in terms of 
its orthogonal components, along some convenient orthonormal basis (We use orthonormal 
wavelets for this purpose). 

Let y B mn ,g e mn and w 9 mn be the discrete w'avelet transform of 9j),g(s t , 0 : ) and w(s,, 0 : ) 
respectively. Because of the orthonormal representation w^ n are independent Gaussian 
random variables. From the linearity property of the wavelet transform it follows that, 

9 6,9 

y mn 9 ran T ^mn 

For notational convenience we have used 6 6 0 in place of 6 : . 

The data available for the detection problem now is {y 9 mn : m € M., n 6 € 

0} where M. — { 1 , 2 , . . . , M} is the set of available distinct difference levels, Af m = 
{0, 1, ... , N m — 1} and N m is the number of samples at each difference level m. We would 
like to determine which of these available detailed coefficients are due to signal and noise 
respectively. The suboptimal test given by equation 2.14 can be used to determine this, 
i.e., if |y^ n | > 7 (m) then it is due to signal, otherwise it is due to noise. The result of this 
hypothesis testing is the detected coefficient set g e mn , which are given by, 


a 9 

9mn 


■I 


Vmn 


^ l^mnl > 7 M 

if \vL 


0 if|»IJ<7(m) (3 ' 1} 

This is precisely the time variant filter which we discussed earlier. The detected coeffi- 
cients g^ n are used to reconstruct the preprocessed signal g(s t ,9j) through inverse wavelet 
transform. This is then used to obtain the noise suppressed image through inverse Radon 
Transform. 
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3.2 Deciding the threshold 

We now turn our attention to the problem of fixing the threshold. In this work, we 
choose a threshold based on energy criterion [11]. As noted earlier, the wavelet transform, 
decomposes a signal into a group of detailed signals, each corresponding to a particular 
frequency band of constant relative bandwidth. The energy at a specific detailed level 
gives the estimate of the signal’s energy in that frequency band. We can estimate the 
information signal energy from the observed signal energy and the noise variance. 

Referring to figure 2.2 and 2.3, we note that with subband implementation scheme the 
detailed signal at the scale 2 m is obtained by passing the original discrete signal through 
m-1 stages of lowpass filters [h(n)] and a single stage of highpass filter [g(n)]. The variance 
of the Gaussian random variable at the output, for N(0,cr 2 ) as input, is given by [8], 

*1 = ^ J2( h ( n )) 2 EteM) 2 (3-2) 

. n J n 

Therefore the noise energy at this difference level is 

= N m ol » 2'“J Noi ■ (3.3) 

Also the observed signal energy ,p Vm , at the level m (scale 2 m ), is given by 

p, m = X>™> ! < 3 - 4 ) 

n=0 

Using equations 3.3 and 3.4 we can estimate the information signal energy, i.e., uncorrupted 
signal energy P g ^ , as 

P g e m = P ym - 2-No-l (3.5) 

This estimate is used to fix the threshold in the equation (3.1). We vary the threshold in 
such a way that the thresholded signal energy P^ = En(?mn) J ) approximates the estimated 
signal energy. 

The proposed technique has been tested on computer simulated images and the results 
are given in chapter 5. 



Chapter 4 

2-D Wavelet Transform 


In this section, we prove that the reconstructed image from the wavelet transform (at a 
fixed scale) of the the projections is the two dimensional wavelet transform of the original 
image at that scale. 

Proof : 

1. The one dimensional wavelet transform of the projection g(s, 0 ), at a particular res- 
olution, a, 

9a{s> d ) = ^j $(*> dt 

can be written as, 

9a(s, 0) = g{s , 0) (~)J 

From the convolution property of the Radon transform (equation 2.6), we have 

G a (x,y) = F(x,y) <g> 2 V a (-x,-y) (4.1) 

where G a (x , y), F(x, y) and 'P a ( — x, —y) are the inverse Radon Transforms of g a (s, 0 ), 
g(s,0) and (1 / y/a)ip(— s/a) respectively; ® 1 ,® 2 denote the one dimensional and two 
dimensional convolution respectively. 

2 . Consider the function, ^ a (— x, — y), which is the inverse Radon transform of — * ), 
i.e., 

'!'„(— x, -y) = - 7 = / / rj>[ — ) p(x cos 0 -f y sin 0 - s) ds d0 

yj CL JO J— 00 V CL / 
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where p(s) is the ramp filter as defined in section 2.1. Since p(s) is symmetric, we 
can write, 

^o(£>y) = —j= J ip ^ p(x cos 6 + y sin 9 — s) ds dd 
With a change of variable, s/a = r, we have, 

d rrt roo 

Va(x,y) = —j= / / ip (r)p(x cos d + y sin 0 — ar) dr dd 

CL J 0 J — oo 

= —7=/ [ ip(r)p a (— cos d + — sin d — r)l dr dd 
y CL Jo J — co . \ (L CL J J 


(4.2) 


We know that p(s) |^|/( 2 x). Now, 


p(as) m \p(s) 

a Ztt a 1 27r cr 

where F.T. and I.F.T. represent Fourier transform and Inverse Fourier transform 
respectively. Hence, 

P (as) = ^p(s) (4.3) 

Using equation 4.2 and 4.3 

i S a (x,y) = - 7 = f f ip(r)p ( — cos 9 + — sin 0 — r\ dr dd 

CL\J CL Jo J— oo \<2 d / 

= _L* 

ay/a Va’ a/ 

where 

/•-r yco 

^(x,y) = / / ip{r)p(x cos 0 + ^ sin 0 — r) dr dd 

Jo J- oo 

Thus we have shown that, 

«.(*,») = ;rs*(M) ( 4 - 4 ) 

fli/a \a a/ 

that is, ^ 0 (a:,y) is the scaled version of the 2D signal $!(x,y). Hence the convolution 
in the equation 4.1 is the convolution with the same signal which is scaled in both x 
and y by an amount a. 


3. Now we show that 'J , (x, 3 /) satisfies the admissibility condition, i.e., 

da < oo 


c , a (2,)’ r Jih^L 1 

J— oo Id] 


(4.5) 
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Now, 


♦to,*) ^ rr^(-) 0K) 


|ap a a |aj 'a' 
where 2dIFT, RT and ldFT represent the 2-D inverse Fourier transform, Radon 
Transform and 1-D Fourier transform respectively. But from Projection slice theorem, 
we know that the 1-D Fourier transform of Radon transform of a function is the same 
as the 2-D Fourier transform of the function. Therefore, 


= fat) 


(4.6) 


From 4.5 and 4.6, we have, 


= 


/«_\2 f °° ^ K ) J n 

(2,r) y-»TT 

2irC^ 


But, from equation 2.8, < oo because the one dimensional wavelet ip(x) is ad- 

missible. Therefore Cy < oo and hence the two dimensional wavelet, ^(x, y) is 
admissible. Here we made use of Fourier Slice Theorem. 


We summarise our result below, 

• Compute the discrete wavelet transform (the detailed coefficients d e mn of the projec- 
tions) and the blurred signal c 9 Mn through subband technique, where m 6 M. and 

U € Njn - 

• Reconstruct the image at a particular resolution from the wavelet transform of the 
projection. 

We give the schematic of the given algorithm in figure 4.1. This property is demon- 
strated using computer calculated projections of simulated images and the results are given 
in chapter 5. 
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e 

d : Detailed signal 

Tfi IV 

(Wavelet transform of 
the projections) 


O 

c : Blurred signal 

Mtv 


B 

F (n,n): Detailed image 
tk 1 x 

e 

B (n,n): Blurred image 
Ml* 


Figure 4.1: Schematic of the given algorithm 






Chapter 5 

Simulation Results 


In this chapter, the results of the chapters 3 and 4 have been verified through simula- 
tion. We first present the simulations results of wavelet preprocessing of noisy projections. 
Then we present the results of the simulations of 2D wavelet transform. We have used 
Daubechies’s compactly supported wavelet of length 16. — 

5.1 Wavelet Preprocessing 

We use two test images to demonstrate the noise suppressing and edge preserving properties 
of the algorithm presented in chapter 3. We add white Gaussian noise to the projections 
of the test image. The thresholding technique(equation 3.1) is used to process the noisy 
projections. These processed projections are then used to reconstruct the image to give 
the noise suppressed image. 

The details of the test image and the simulation results are given below. 


• a. 


< 0.5 and |y| < 0.5 
otherwise 


/(*.»>-{ j “ w 

The function is defined over the square —1 < x,y < 1. The Radon transform, g(s,0), 
of this function, for — \/2 < s < -y/2 and 0 < 8 < ir/4 is given by, 


9(- s,6) 


( j+o.5cosg+o.5sin? —0.5 cos 9 — 0.5 sin 6 < s < —0.5 cos & + 0.5 sin 0 

1 sin 9 COS 9 

l 

cos 6 . 

— j-f-o 5 cqs$+q s sm § _ q . 5 cos 0 — 0.5 sin 8 < s < 0.5 cos 0 + 0.5 sin 0 

sin 9 cos 9 


—0.5 cos 0 + 0.5 sin 0 < s < 0.5 cos 0 — 0.5 sin 8 


Symmetry of the figure is used to obtain the value of g(s, 0) for 7r/4 < 0 < it. 
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We have used samples of g(s,6) with 127 samples in variable s (number of rays per 
projection) and 40 samples in variable 9 (number of projections) to reconstruct the 
test image on a grid of 51 x 51. Reconstruction from noiseless projections, g(s,9), is 
shown in figure 5.3. We add white Gaussian noise with zero mean and variance 0.0004. 
Reconstruction from these noisy projections is shown in figure figure 5.4. The noisy 
projections are then processed using the thresholding technique(equation 3.1). The 
reconstruction from the processed projections is shown in figure 5.5. Reconstruction 
from noisy projection was also done using the space invariant filter(equation 2.4) 
technique. To compare the reconstruction using the thresholding technique with the 
space invariant filtering technique, we plot the cross section through each image along 
the line x=0 in figure 5.7. From figure 5.7, we observe that the reconstructed image 
after preprocessing closely approximates the noiseless reconstructed image at the 
edges. However, we observe from figure 5.8, this is not the case for the reconstructed 
image using space invariant filtering technique. We also note from figure 5.7 that 
the noise performance of the thresholding technique is as good as the space invariant 
filtering technique. 

• b. Shepp-Logan’s Head Phantom : 

We now repeat the above process on Shepp-Logan phantom [13]. This phantom, 
which is frequently used in medical imaging, is supposed to be a cross sectional view 
of the human head with denser areas indicating tumors and the less dense areas in- 
dicating the spinal fluid [9]. Figure 5.1 shows this phantom. We have used samples 
of the Radon transform g(s, 9) with 255 samples in the variable s and 100 samples 
in variable 9 to reconstruct the test phantom on a grid of 99 x 99. Reconstruction 
from noiseless projection is shown in figure 5.9. We then add white Gaussian noise 
with zero mean and variance 1 x 10 -6 . Reconstruction from these noisy projections 
is shown in figure 5.10. The noisy projections are processed using the thresholding 
technique. Image reconstructed from these processed projections is shown in fig- 
ure 5.11. Figures 5.13 to 5.14 show the cross sections of these reconstructed images 
along the lines x=0 and y=0. 



Chp.5 Simulation Results 


21 


5.2 2D wavelet transform 

The results of chapter 4 are demonstrated on a test image. The test image used in this 
case is shown in figure 5.2. Reconstruction is done on a grid of 79 x 79 from 50 projections 
with 255 rays per projections. The reconstruction from noiseless projections is shown m 
figure 5.15. We consider the original discrete projections to be at resolution level 2°. Using 
subband technique (refer figure 2.2), detailed signals at resolution levels 2 -1 and 2 2 (i.e. 
the wavelet transform at the scale a=2 and 2 2 ) of the projections are calculated along 
with the blurred signal at 2~ 2 . Reconstructed images from the detailed signals and the 
blurred signal of the projections at these resolution levels are shown in figures 5-16-5.18. 
We note that the image reconstructed from the detailed signal give the local high frequency 
information i.e., the edge information, of the original image. 

We add white Gaussian noise to these projections with variance 0.0001- After thresh- 
olding is done in the wavelet transform domain, instead of taking inverse wavelet transform 
to get the processed projections, we reconstruct the images as mentioned before. This re- 
sults in the noise suppressed 2 dimensional wavelet transform of the original image from 
noisy projections. The image thus reconstructed from the noise supressed detailed signal, 
at resolution 2 _1 , of the projections is shown in the figure 5.19. The wavelet transform at 
this scale represent the high frequency part of the projections. Thus here noise dominates 
the signal. In figure 5.19, we observe that the the thresholding technique has essentially 
removed this noise. 
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Figure 5.1: Shepp-Logan Phantom, 
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0.0000 
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0.6624 
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d 
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0.1600 
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g 

H 


0.0460 

0.0460 
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0.01 
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-0.08 

-0.6050 

0.0460 
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0.00 

0.01 
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-0.6050 

0.0230 

0.0230 

0.00 

0.01 

| j 

0.00 

-0.6060 

0.0230 

0.0460 

0.00 

0.01 


Tkble 1: The parameters of the ellipses in Shepp- Logan's head phantom. (r,y), 
roDresents the centre of the ellipse. 
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Figure 5.2: Test Image (Elipse and Circle) 
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Figure 5.7: Cross section along x=0 of image reconstructed from wavelet preprocessed pro 
jections 
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Figure 5.8: Cross section along i=0 of image reconst, acted from ..pact invariant Jilted 
projections 
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Figure 5.9: Reconstruction from Noiseless projections 



Figure 5.10: Reconstruction from 
Noisy projections (Noise Variance 

i x io" 6 ; 


Figure 5.11: Recon- 

struction from Preprocessed projec- 
tions (Noise Variance 1 x 10~ 6 J 
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Figure 5.12: Cross section of phantom reconstructed from Noisy projections (Noise Van 
ance 1 x 10~ 6 j 
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Figure 5.13: Cross section of phantom reconstructed from Preprocessed projections (Noi. 
Variance 1 x 10~ 6 J 



Chp ) ^iii hildt ion Results 


.{U 



Figure 5.14: Cross section of phantom reconstructed from space invariant filtered projec- 
tions (Noise Variance 1 x 10 -6 ) 
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Figure 5.17: Detailed image at reso- 
lution 2 _1 


Figure 5.18: Detailed image at reso 
lution 2 -2 
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Figure 5.19: Noise suppressed Detailed image at resolution 2 1 (Noise variance 4 x 10 4 ) 



Chapter 6 

Scope for Future Work 


6.1 Reconstruction of 2D image from its Wavelet 
transform 

In chapter 4, we outlined a method to obtain the two dimensional wavelet transform of 
an image from its projections. Here we suggest a method to synthesise the image from its 
available wavelet transform. 

From equation 2.12 we have, 

h(n ) = -^= J ^(s/2)?!i(s — n) ds 
= ^~ s ) 

where n G Z. Let h c (s ) = (l/\/2)(f>(s/2) <gq <f>(— s), i.e., h(n) is the discrete signal obtained 
by sampling h c (s) at 5 = n;n € Z. Let H c (x,y) represent the inverse Radon transform 
(IRT) of h c (s). From convolution property (equation 2.6) and the scaling property [1], 

where $(x, y) is the IRT of ^(s). Similarly from equation 2.13, we can write, 

G c (x, y) = ® 2 _y) 

where G c (x,y) is the IRT of y c { s ) (which can be sampled at s = n\n € Z to give g(n)) 
and tf(z,y) is the IRT of 1>(s) (equation 4.4). Let H(n u n 2 ) and G(n u n 2 ) represent the 
discrete forms of H c (x,y) and G c {x,y) respectively obtained by sampling at x = n u y - 
n 2 ; , n 2 G Z. 
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Let F m (ni,n 2 ) and B m (ni,ra 2 ) represent the detailed image and the blurred image at 
the resolution 2 obtained by the method suggested in chapter 4. We propose that the 
blurred image i?( m -i)(fti, n 2 ) at the next higher resolution viz. 2 _ h n ~ 1 ) can be obtained as 
follows, 

• Upsample both F m (n u n 2 ) and B m {n u n 2 ) by a factor of 2 in both the directions n 2 
and n 2 . 

• Filter the upsampled blurred image with f/(n 2 , n 2 ) and the upsampled detailed image 
with G'(n 1 ,n 2 ). 

• Add the filtered images. 

We justify our claim by repeatedly invoking the linearity, convolution and scaling prop- 
erties of Radon transform. Let c 8 n and d 8 ln represent the blurred and detailed signals at 
the resolution 2 _1 of the projections of the image Ro(fti,rc 2 ). These signals are obtained 
by filtering the projections using h(n) and g(n) respectively and down sampling them by 
2. Figure (6.1) shows how to reconstruct the image from and 

• Making use of linearity property of Radon transform, we can shift the Inverse Radon 
transform (IRT) block before the adder, noting that the addition will now be two 
dimensional. 

• Making use of convolution property, we can shift the IRT block before the filtering 
stage. Now the filtering will be two dimensional and it is the filtering of the re- 
constructed blurred (or detailed) image with the IRT of the function whose discrete 
representation is h(n) (or g(n) as the case may be). 

• Making use of scaling property, we can shift the IRT before the downsampling stage, 
noting that the down sampling will now be done in two dimensional and will be in 
both direction and n 2 . 

The block diagram after the above modification is as shown in figure 6.2. Since the blurred 
signal cj n can be obtained from c 2n and applying the same 3 steps, we can show that 
the IRT of the blurred signal c? n and hence the image, B 0 (ni,n 2 ), can be obtained from the 
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h(n) : Lowpass filter c : blurred siqnal 

B (n } n) : Blurred image 

g(n) : Highpass filter d : detail signal 
( 1-D filters) (1-D signals) 

Figure 6.1: Reconstruction of an image from the detailed and blurred signals of its projec- 
tions 





H (n,n) : S-D Lowpass filter 


( refer text for details) 


G (n,n) : S-D Highpass filter 

Figure 6.2: Modified block diagram 


IRT of c 2n and d e 2n . In the same way, we can obtain the IRT of the blurred signal c* m _ 1)n 
which is the blurred image n 2 ) at the resolution can be obtained from 

the IRT of the blurred signal c 6 mn , (the blurred image 5 m (nj,n 2 )) and the IRT of the 
detailed signal d 9 mn (the detailed image F m (ni,n 2 )) at resolution 2~ m . 

Thus we have outlined a method to hierarchially reconstruct an image from its two 


dimensional wavelet transform. 
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6.2 Reconstruction of localised area of interest 

An important application of multiresolution reconstruction of the tomographic images, is in 
reconstructing only the localised region of interest at high resolution with the sorrounding 
area being reconstructed at a lower resolution to provide a frame of reference. Such a 
technique may be useful in medical imaging, where a physician may be interested in the 
cross section of a small region [18]. In [17], Olson et.al., present an algorithm wherein, 

• projections at certain angles are exposed only over the localise area of interest; 

• one dimensional wavelet transform is applied to each projection to find all of the 1-D 
wavelet coefficients; 

• the low resolution wavelet coefficients for the reduced-exposure projections are re- 
placed with coefficients found by angularly interpolating between the low-resolution 
wavelet coefficients of the full exposure projections; 

• one dimensional inverse wavelet transform is then used to create a complete set of 
modified projections; 

• these modified projections are then used in a standard reconstruction algorithm to 
arrive at a reconstruction which has high resolution in the region of interest. 

The key to their algorithm is that the space-frequency localisation property of compactly 
supported wavelets (e.g., Daubechies [5]) is essentially preserved after applying the ramp 
filter used in FBP algorithm [17, 18]. The major drawbacks in their algorithm are listed 
below. 

1 . The angular interpolation involved introduces error in the reconstruction [18]. Though, 
it may not affect the reconstruction at the region of interest, the reconstruction have 
severe artefacts outside the region of interest. 

2. Their algorithm is prone to aliasing error when the region of interest is off centered. 
The process of centering the projections before interpolation and reconstruction in- 
troduces this aliasing error [17]. 
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Reconstructed image 
with high resolution 
tn the region of interest 



N : Total number of projections 

Half the number of projections (Reduced exposure Projections) 

are exposed only over the region of interest , t 

Figure 6.3: Schematic of the proposed algorithm for 2 levels. The reduced exposure projec- 
tions are interlaced between the full-exposure projections. 

3. They have not mentioned in their paper how multiple non concentric regions would 
be handled. 

We suggest a modification in their algorithm to overcome these drawbacks. We first 
note that the main reason for all these problems is the angular interpolation involved. We 
follow the same radiation exposure pattern suggested by Olson et.al. in [17]. We eliminate 
the angular interpolation by reconstructing the images at different resolution level from 
the wavelet transform of the projection using the technique presented in chapter 4. We 
have proved in chapter 4 that these reconstructed images are the two dimensional wavelet 
transform of required image. Then we hierarchially reconstruct the final image from their 
wavelet transform through the method presented in section 6.1. The block diagram of 
the suggested method for 2 levels is shown in figure 6.3. For a schematic of the radiation 
exposure pattern refer [17, 18]. The verification of these two suggestions through simulation 
and if the method is successful, then providing the proper mathematical foundation for it, 
can be taken up in the future. 








Chapter 7 
Conclusion 


In this thesis work, we had studied some applications of wavelet transform in the area of 
image reconstruction from projections. An algorithm was proposed to suppress the effect 
of noise in the reconstructed image from noisy projections using wavelet transform based 
on spatially variant filtering technique. We demonstrated through simulations, the edge 
preserving property and the noise suppressing capacity of the proposed technique. We 
presented a qualitative comparison of our technique with the existing technique on space 
invariant filtering techniques. We also proved that the images reconstructed from the one 
dimensional wavelet transform of the projections at a given scale are the two dimensional 
wavelet transform images of the original image. The two dimensional analysing wavelet is 
the inverse Radon transform of the 1-D wavelet. The algorithm was implemented to obtain 
the two dimensional wavelet transform images. 

As an extension to the present work, two suggestions were made in chapter 6 for the 
future work. One deals with reconstructing the image from its wavelet transform obtained 
by the above method while the other deals with the reconstruction of specific region of 
interest from essentially localized projections. 
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Appendix A 


We attempted to quantify the error in the reconstruction from noisy projections with the 
preprocessing technique. The results are presented here. Two error measures are used [3]. 

MSE= _t ( m ’ n )] 2 ) • (A- 1 ) 

\ m n } 


NMAD = 1*2. !/(”*,") -*(™,n)| 
Em Ej*(m,n)| 

In the above equations, 


(A.2) 


MSE 
NMAD 
/(m, n) 
t{m , n) 
t 


Mean Square Error 
Normalised Mean Absolute distance 
Reconstructed image 
Original image 

Average value of the original image 


The simple square image which is given in the simulation chapter is used. The MSE and 
NMAD as a function of the variance of the noise added to the projections is plotted in 
figure A.l and A.2 for the images, 


• reconstructed from noisy projections. 


• reconstructed using hamming filter technique. 

• reconstructed using the proposed preprocessing technique. 

The variance of the noise in the projections is varied from 1 x 10 -4 to 36 x 10 -4 The error 
in the image reconstructed from the hamming filter technique is large due to the large 
deviations near the edges in the image. 
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{ Haondng filter 



Figure A.l: MSE vs Noise variance 
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Figure A.2: Normalised Mean Absolute Distance vs Noise variance 



